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ABSTRACT 

We make a revision of the stability criteria for equatorial circular orbits, obtained from the epicyclic 
approximation, which is widely used in Newtonian models for axisymmetric galaxies. We find that, 
for the case of thin disk models, the criterion for vertical stability must be reformulated, due to the 
discontinuity in the gravitational field. We show that, for a model characterized by a surface mass 
density S, the necessary and sufficient condition to have vertically stable circular orbits is that £ > 0. 
On the other hand, the criterion for radial stability is the same as in thick diks, i.e. that the (radial) 
epicyclic frequency squared is positive. As an application, we present finite thin disk models for nine 
galaxies, as superpositions of members of the Morgan & Morgan family (in Newtonian version), which 
can be considered as stable configurations in a first approximation. Also, as an additional product 
of this study, we show that any galactic model with a thin disk component admits a wide variety of 
integrable disk-crossing orbits, which are determined approximately by a third integral of motion of 
the form ZE 1 / 3 , where Z is the z-amplitude of the motion. 
Subject headings: stellar dynamics - galaxies: kinematics and dynamics. 

1. INTRODUCTION 

It is usually accepted that many galaxies in the universe are nearly axisymmetric, with a mass distribution formed 
by several components: a thin disk, a central bulge and a surrounding halo. In con sequence, there i s a number of mas s 
models incorporating one or all of these features, depending on the particular case (jFreemanI ([l97f)h : iKentl (fl98l Il987h : 
iSofue. Honma fc Omodakal (|2009f >). For example, there exists a number of galax ies in the Ursa Major cluster which can 
be mo deled, at the scale of the optical radius, with only the thin disk comp onent ([Gonzalez. Plata-Plata fc Ramos-Carol 
(2010)), suggesting that they obey the so-called maximum disk hypothesis (|Binnev fc Tremaind (|2008l )L If we decide to 
continue using Newtonian theory beyond the optical radius, presumably we have to introduce a dark halo, but the disk 
component still provides a significant contribution, taking into account that the main part of the stellar population is 
located there . For this reason, thin disk models have been an issue of interest in galactic dynamics (see for e xample 
iHvmterl (11963ft: iMorgan fe Morgan! (|1969l ): lGonzalez fe Reinal (|2006f ): IPedraza. Ramos Caro fe Gonzalez! (|2008[ ) as well 
as iBinnev fe Tremaind ()2008l ) and references therein) . 

Thin disks have also been used to model self-gravitating rings ()Lemos fc Letelierl (fl99l : ILetelierl (l200l L with 
applications in other branches of ast rophysics. These m o dels c a n be used to describe accretion disks, st e llar s truc- 
tures with central black holes (see iSemerak fc Sukoval (|2010( k lLora-Claviio. Ospina-Henao fc Pedrazal (l2010f)) o r 
planetary rings. For example, the last issue was partially encompassed in fRamos-Caro. Pedraza fc Letelierl |2011| ). 
where a study of the linear stability of the monopole-ring system was performed by using superpositions of Mor- 
gan fc Morgan disks. Similar studies w e re conducted by the authors in the Newtonian realm of galactic dynamics 
(|Ramos-Caro. Lopes-Suspes fc Gonz alez (200~8() : IPedraza. Ramos-Caro fc Gonzalez! ()2008D ). 

A fundamental step in the formulation of galactic models is the stability analysis. In fact, the stability analysis 
could suggest sometimes the introduction of new features in a given model in order to, p r esuma bly, obtain a more 
realistic representation. Consider for example the study conducted by lOstriker fc Peeblesl (|1973| ) where it is shown 
that a flattened system of self-gravitating particles, initially supported against gravity by rotation, does not maintain 
its discoidal form in the course of time. They suggested that the introduction of a spherical halo with mass of the order 
of the disk mass (or greater) improves substantially the stability of the disk, as it was corroborated by simulations. 
But, on the other hand, consider also the counter argument made by iKalnali (|1987l) . in which it is suggested (i) 
that the stability problems can be overcome by improving the features of the the inner parts of galactic models (for 
example, by considering hot centers or small bulges) and (ii) that a halo with a scale length larger than the disk and 
more massive than it, does not contribute significantly to the stability. This discussion has its roots in the fact that 
internal kinematics of self-gravitating disks determi nes the stability of the syste m: cool disks, which are mainly formed 
by stars in circular motion (the ones considered by lOstriker fc Peeblesl (11973|)). requ i re a p rominent halo in order to 
avoid instabilities, contrary to the case of hot disks (the ones considered bv lKalnais! (1987)) which can be supported 
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by the random motions of the stars (for a recent review of disk instabilities see ISellwoodl (golil)). A detailed knowledge 
of the orbital features associated with a given galactic model provides the basis for the stability analysis. Usually this 
knowledge can be achieved once we have at hand the distribution function (DF) of the model. 

In general, the obtention of the DF associated to a particular model is not an easy task (except for those models which 
are defined by an analytical DF) and, in consequence, the corresponding stability analysis is far from trivial. However 
it is possible to perfo r m a fi rst test of stability without knowing the DF, by using the so-called epicyclic approximation 
(|Binnev Sz Tremain c (2008)), i.e. by performing a linear stability analysis of fixed points of the effective potential. 
Once it is verified that the model is linearly stable, it deserves to perform more conclusive stability tests based on 
statistical mechanics. 

The linear stability analysis can be thought as a first test to evaluate how realistic any particular model is. In disk 
galaxies many stars are on nearly circular orbits, so we have to demand, as a basic requirement, that any galactic 
model must be characterized by allowing the existence of stable circular motion, especially in regions where the stellar 
population is maximum, i.e. the equatorial plane. The epicyclic approximation provides a formalism to study motion 
in the equatorial plane and leads to the establishment of simple stability criteria. However, when one deals with models 
incorporating a razor-thin disk, this method needs to be reformulated in view of the mathematical features introduced 
by the thin disk component. Next, we will briefly illustrate this problem. 

Consider an axisymmetric potential-density pair (APDP) with mass distribution p(R, z) and gravitational potential 
$(i?, z), wher e R and z are the usual c ylindrical coordinates. The motion of test particles can be described by the 
Hamiltonian ()Binnev fc Tremaind (|2008D ) 

p2 , p2 

H=-?^ + <S> eff (R,z), (1) 
where Pr = R, P z = z and & e ff is the so-called effective potential, defined as 

* eff (R,z)=*(R,z) + -±p. (2) 

Here, I represents the z-component of the angular momentum, often called azimuthal angular momentum, which is a 
first integral of motion. The resulting equations of motion can be written as 

<*> 

In particular, the equatorial circular orbits (i.e., belonging to the plane z = 0) correspond to the minimum of $ e //, 
obtained by setting d& e ff/dR = and z = 0, which lead to the relation 



<9$ 
dR 



I 2 

-S»=0- (5) 



The value of the i?-coordinate that is solution of the above equation is the radius of the circular orbit with angular 
momentum I (in this case R is called the guiding- center radius, which we will denote as R ). If this orbit (or any other 
in the equatorial plane) suffers a small perturbation, one would expect that the resulting motion is not very different 
from the original one (say, a nearly circular orbit), in order to guarantee the stability of the entire configuration. 

Strictly speaking, a nearly circular orbit is defined as a trajectory with coordinates (R, z) very close to (i? o ,0), so 
we can expand <& e ff near this minimum and neglect cubic and higher terms in the expansion, that is, 



where k and v are defined as 



K 2 d^ eff 



dR 2 



v 2 _ d 2 ^ eff 



(flo,0) 



(7) 

(Ra,0) 



By introducing (|6]) in the equations of motion ([3])-(|4]), it can be seen that R — R and z evolve like the displacements 
of two harmonic oscillators with frequencies k (epicyclic frequency) and v (vertical frequency), respectively, once it is 
guara nteed that k 2 > and v 2 > 0. In this case the corresponding circular orbit is said to be stable (|Binnev fc Tremaind 
(2008)). So we say that, in a first approximation, the conditions for stability of the self-gravitating structure located at 
the equatorial plane (which it is assumed to be composed principally by particles describing nearly circular motions) 
is that the epicyclic and vertical frequencies squared are both positive. 

Note that, whenever one deals with razor-thin disks, it is not possible to define a Taylor expansion of the ef- 
fective potential around the point (R o ,0). This procedure whic h leads to the analysis of vertical and epicyclic 
frequencies around the circular orbit (Bin nev fc Tremaind (j2008) ) , is only valid for potentials that are smooth 
or that have at least continuous second derivatives. This means that the analysis of vertical frequencies of thin 
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disk potentials, as did in iPedraza. Ramos-Caro fc Gonzalez! (120081); IRamos-Caro, Lopes-Suspes fe Gonzaled (|2008|) 



iGonzalez. Plata-Plata fe Ramos-Carol ( 20101 ): IRamos-Caro Pedraza fe Letelierl ( 20111 ). is not a reliable indicator of 
the vertical stability of the corresponding circular orbits. The analysis must take into account the discontinuity in the 
partial derivative of the potential with respect to z due to the surface distribution of matter in the equatorial plane. 

In the following sections we will show that, when considering a razor-thin disk, it can be constructed a vertical 
stability criterion in terms of the first derivative of the potential, whereas the criterion of the epicyclic frequency 
remains unchanged (sees. [2] and [3]). Also we will present some consequences of this approach in the description of 
disk-crossing orbits. In particular we will show that many of them can be described by an approximate third integral 
of motion depending on the vertical amplitude and the surface mass density (sec. |4}. Finally we briefly address these 
issues in the realm of modified theories of gravity (sec. EJ. 

2. GALACTIC MODELS VIA THIN DISKS 

As it was pointed out in the Introduction, many galaxies are modeled as an axisymmetric thin disk surrounded by 
an axisymmetric 3D smooth density distribution, which is symmetric with respect to the equatorial plane. The total 
density distribution can be written as 

p(R,z) = Z(R)6(z)+p s (R,z), (8) 

where S is the Dirac delta, S(i?) is the surface density distribution of the thin disk and p s (R,z) describes the sur- 
rounding matter. The gravitational potential of the system is 



$(R,z) = $ E (i?,z) + $ s (i?,z), 
where $ s and $e are the contributions due to p s and S, respectively, in such a way that 



(9) 



V 2 $ s = 47rGp s , V 2 $ E = AttGH{R)5{z) 



The symmetry of p s with respect to the equatorial plane implies a reflection symmetry of the gravitational potential 
with respect to the equatorial plane, i.e. $(R, z) — <t(R, —z). This reflection symmetry implies that the z-dependence 
of $ actually is solely on \z\, so we can write 



<9$ _ d<I> d\z\ _ <9$ 
~d~z~ d\z]^~ WY 



(10) 



where O is the Heaviside step- function. Thus, the potential $ in these models is a smooth function of \z\ with a 
discontinuity in the z-derivative at the plane z = 0. 

It is suggestive to note that there is a close connection between d$/d\z\ and E. To see this, first note that a direct 
consequence of (flOl) is 



9$ 

~dz~ 



o+ 



<9$ 

dz 



Now, from Gauss's theorem we have 



and, according to (fTTj). it follows that 



47rGS(i?) 



9$ 

dz 



= 2 



o+ 



<9$ 

W\ 



(11) 



z=0 



dz 



1 d<S> 



w=^cm^ (12) 

which is equivalent to the expression appearing in Gonzal ez fe Reinal (|2006D and references therein. This fact will be 
relevant in the construction of a stability criterion for circular orbits under the action of vertical perturbations. 

3. VERTICAL STABILITY OF CIRCULAR ORBITS 

We pointed out in the Introduction that the stability analysis of equatorial orbits must take into account the 
discontinuity in the partial derivative of the potential with respect to z, due to the surface distribution of matter in 
the equatorial plane. The aim of this section is to analyze the behavior of the z-coordinate of a vertically perturbed 
orbit in the vicinity of the discoidal distribution. 

Consider an equatorial circular orbit of radius R under the action of a small vertical perturbationQ. This perturbation 
can be seen as an instantaneous vertical increase, v 0z , in the velocity of the particle, which does not affect the value of 
/. In order to study the evolution of the perturbation in the course of time, we have to use the z-equation of motion. 
From eqs. (0| and (ITUl) . we have 

z = -[ 2 e(z)-l]|i (13) 



from which we can establish the following statements: 



1 Here, the term "small" means small enough to neglect variations in the projection of the orbit on the z = plane. 



4 



1. Suppose that the particle hits the disk from below. If 



d<S> 



> then 



< z 



therefore its vertical acceleration decreases and the particle tends to come back to the disk. 



2. Suppose that the particle hits the disk from above. If 



> then z 



> z 



therefore its vertical acceleration increases and the particle also tends to come back to the disk. 



< 0. 



z=0 



3. Consider one of the two above situations but with 

W\ 

In both cases, the change in acceleration due to the discontinuity in the z-derivative of the potential makes the 
particle to gain velocity and, intuitively, move out from the disk. 

From the above considerations, it is natural to establish that a necessary condition for vertical stability is given by 
the relation 

> 0. (14) 

z=0 



a\z\ 



But the vertical stability is guaranteed once we verify that any particle that crosses the disk, with initial vertical 
velocity Vq z (just after crossing), will come back to the disk. We show in the appendix that, for vq z sufficiently small, 
the perturbed trajectory will oscillate around the original one for large time. Then we can consider that condition 
(fl~4"|) is also sufficient to guarantee vertical stability. If we start with a small vo z , the particle will always cross the 
disk with a velocity whose vertical component has modulus \vq z \. This is a consequence of: (i) the conservation of 
mechanical energy; (ii) the assumption that the incoordinate does not change in the process and (iii) the fact that the 
discontinuity is only in acceleration. So, just after crossing the disk the particle will also have a velocity with vertical 
component of magnitude Vq z , which means that the motion after crossing the disk will have the same behavior as the 
motion before crossing it. Therefore, the perturbed orbit remains oscillatory around the original one for sufficiently 
small initial vertical velocity vq z ■ 

According to (fl"2]) . condition (| 14[) can be written in terms of the surface mass density. Then we can state that in 
a Newtonian thin disk model, a necessary and sufficient condition for a circular orbit of radius R to be stable under 
small vertical perturbations is 

E(R) > 0. (15) 

This condition has indeed the same status as the vertical stability condition for smooth potent ials (that is, replacing £ 

by v 2 in the above statement, as was employed in lGonzalez. Plata-Plata fc Ramos-Carol (f2010l ); lRamos-Caro. Pedraza fc Letelierl 
pOlM in the sense that both conditions imply vertical stability of the circular orbit if we neglect variations in the 
i?-coordinate of the perturbed orbit (see appendix) . 

In particular, if there are regions where E(i?) = 0, i.e. without a thin disk component, then in these regions we must 
apply the criterion of vertical frequencies to analyze the vertical stability of the corresponding circular orbits. 

For the case in which the perturbation introduces small variations in the ^-coordinate, we must know something 
about the radial dependence of the effective potential to ensure stability of the perturbed orbit. It turns out that a 
sufficient condition for the i?-variations in the perturbed orbit to be negligible is (see appendix) 



k 2 (R) > 0. 



(16) 



We also point out that the simultaneous conditions l|15 p -(|16 p imply (Liapunov) stability of the circular orbit under small 
perturbations in an arbitrary direction of the meridional plane (by neglecting variations in I due to these perturbations). 
We will address this issue in detail in Sec. [4j 

By way of verification of condition 1|15|) . in the following subsection we show that, for a sufficiently small vertical 
perturbation, it is possible to compute characteristic periods and amplitudes of the oscillations around the discoidal 
plane. 

3.1. Characteristic period and amplitude of the oscillations 

Without loss of generality, consider a circular orbit of radius R which suffers a vertical perturbation in the direction 
z > at time t = and neglect its radial variations. In this moment, the test particle has an initial vertical velocity 
vq z > and starts to rise, but it is attracted to the equatorial plane by gravity. The equation of motion for the vertical 
coordinate, restricted to the region z > 0, is 

r t 

z(t) - v 0z 



Suppose that after a time interval At the particle returns to the disk. Then we have 

(-Vox) ~ V0z 



s to the disk. Then 
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For sufficiently small vq z we can approximate the integrand by a constant, obtaining 

2v Oz *(At)^-(R,0+). 

oz 

We can perform an analogous procedure for the region z < and the result is the above equation but with a minus 
sign. Then, we can extend the above relation for any real value of z: 

2v 0z * A* 7— (R,0). (17) 
d\z\ 

Since the characteristic period of one oscillation is T = 2At, eq. (|12p gives us 

T =^m-y < 18 > 

The amplitude of the oscillation is given by the maximum of the parabola described by the equation 

Z(t) S3 V 0z - t——r 

d\z\ 

and, since z(0) = 0, we have 

t 2 d<$> 

Z(t) S3 V 0z t 



z=0 



2 8\z\ 



z=0 



v 2 = V 2 $ 



The maximum of z, which we will denote as z max , occurs at t = T/4, and is given by 

v 2 

Zmax = 4nGZ(RY (19) 

As expected, the amplitude of oscillations is inversely proportional to the surface mass density (but not to the volumetric 
mass density, od the surrounding matter, if present), a natural consequence of the attractive gravitational force. These 
estimates give the order of magnitude of the parameters. More rigorous and accurate estimates are discussed in the 
appendix. 

3.2. Thin disks as a limiting case of a high density thick disk 

The stability condition obtained makes sense physically. For a thick disk, the vertical stability of a circular orbit is 
given by the sign of the square of the vertical frequency, 

1 dv 2 

( 2 °) 

z=0 R dR 

In general, the density distribution is assumed to be very high around the plane z = and to fall off quickly with z. 
If we consider that the rotation curves are bounded, the last term in the above equation is also bounded (remember 
that v 2 — Rd$/dR), and the sharpening of the density distribution makes the first term on the right-hand side of (|20[) 
increase. Eventually, it will increase an amount such that v 2 is positive (this may depend on R). In the limiting case 
when the (positive) density distribution becomes very close to S(z), v 2 will be positive everywhere on the disk. 

3.3. Finite Thin Disk Models for Galaxies 

Some authors have demonstrated that it is possible to obtain galactic models with the thin disk component only 
(without invoking a halo) , model ing with high p r ecision the r o tation curves in the optical range and in acco rdance with 
reasonable mass density profiles (jKalnaisI d 19831 ) ; lKentl (|1986[ ); IGonzalez. Plata-Plata fc Ramos-Carol ([201 0)). Galaxies 
that can be described by these models, which usually have high surface brightness (HSB), ar e said to obey the so- 
called maximum disk hypothesis (|Binnev fc Trerna inc (2008)). This property was confirmed bv lPersic. Salucci fc Stell 
(|1996f) , where they found by means of mo del- independent methods that the amount of dark matter in HSB galaxies is 
negligible insid e the optical radius. 

In particular, iGonzalez. Plata-Plata fc Ramos-Carol ([2010) show some examples in the Ursa Major cluster, by using 
the Hunter method 0. As it was pointed out, one of the conclusions stated by Gonzalez et al. is that the models 
obtained are vertically unstable, due to the fact that they considered the stability criterion of the quadratic vertical 
frequency. But this is not true, in light of the stability criterion constructed here and taking into account that all of 
the mass density profiles obtained in the aforementioned reference are positive. In fact, we remark that the formalism 
showed by Gonzalez et al. is a powerful method to obtain stable maximum diks describing the optical region of many 
HSB galaxies. 



2 This is the same formalism to obtain the Gener alized Kalnajs Disk s (Gonzalez & Rcina (2006)), which are the Newtonian version of 
the general relativistic Morgan & Morgan solutions (Morgan & Morgan (1969)). 
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TABLE 1 Parameters for the nine spiral galaxies chose n: a is the radius of the optical disk, m is the nu mber of ex pan sion 
constants used to obtain the best fit (i.e. the A^n of eq. ((2lJ), C^n are obtained from the fits a ccor ding to (I2.'il and 1241) . The 
estimates for the total mass of the maximum disk are computed from (125 D . 

In order to illustrate the above statements, we present additional models for nine spiral galaxies: NGC 55, 
1417, 3495, 3672, 3691, 4062, 4605, 5585 a nd 5907. It can be assumed that all of them obey the rotation law 
(jGonzalez. Plata-Plata fc Ramos-Carol (j2010f )) 



(21) 



where a is the radius of the stellar disk, i.e. the optical radius, and A 2n , along with the integer m, are constants to 
be determined by fitting the observational data corresponding to the circular velocity v c . The surface mass density of 
the thin disk is given by the relation 



E ( fl ) = 7T^y) C 2n(2n+l)g 2n+1 (0)P 2n (77), 
znaGrj 



(22) 



where r\ — yjl — R/a, P n are the Legendre polynomials, q n (£.) = i n+1 0nK), Qn being the Legendre Functions of 
second kind, and C2n are determined from Am through the relation 



C*2n — 



An + 1 



4n(2n + l) 92 „(0) 



for n and 



m „i 

Y^Mk / x(l~x 2 ) k Pm(x)dx, 
fc=i 



C = ^(-1)" +1 C 2 „, 



(23) 



(24) 



from which the total mass of the stellar disk can be determined: 

M = aC a /G. (25) 

In figure [T] we show the fits to observational data, corresponding to rotation curves from ISofue etail (fT999K The 
results of the fitting, i.e. the values of m and A 2 „, determine the constants C 2 „ (see table [lj and the mass density 
profiles (eq. (|22|) '). which are shown in figure[2l In addition, we can obtain the quadratic epicyclic frequency (fig. [3]) 
from the relation 



K 2 (i?) = ^2(n+l)A 2 „(i?/a) 2 "- 2 , 



(26) 



which is provided by Gojnzalez^, Plata-Pl ata fc Ramos-Carol (|2010l ) . We can see that in all of these cases it is possible 
to obtain models whose circular orbits are all stable. 
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Fig. 1. — Circular vel ocity, V c (km s ) , vs. R/a (optical disk region) for a sample of n ine spiral galaxies. The points correspond to 
observational data from Sofuc ct al. (1999) and the full line is the best fit using relation 121> . 

4. INTEGRABILITY OF MOTION NEAR A STABLE CIRCULAR ORBIT 

We remarked in the Introduction that motion is integrable near a stable circular orbit. This fact follows from the sep- 
arability of the Hamilton- Jacobi equation near (R , 0). There a re many evide n ces of this general behavior for thin disks 
in the l iterature dGonzalez. Plata-Plata fc Ramos-Carol f2010D;lHunterl (120051) ;lRamos-Caro. Lopes-Suspes fc Gonzalezl 
(|2008[ k lSaa fc Venegerolesl (|1999f ): lRamos-Caro. Pedraza fc Letelierl (|201lD : lPedraza. Ramos-Caro fc Gonzalez! (|200cf) )T 
where it is found numerically, by means of Poincare sections, that motion is integrable around what appears to be 
a stable point of the effective potential, corresponding to a stable circular orbit in the equatorial plane. These evi- 
dences also show that the integrable region goes well beyond the neighborhood of the stable circular orbit, where the 
approximation of a separable potential not valid. 

As in the case of a smooth density distribution, for density profiles of the form ([8]) motion near a stable circular 
orbit of radius R in the plane z = is nearly integrable. This follows from the separability of the effective potential 
near the stable point {R D , 0). In fact, up to first order in z, we have 

$ eff (R, z) « $ e// (i?, 0) + -q$-(R, 0)\z\ 

= <$> e ff(R,0) + 2<n:GE(R)\z\. (27) 

Since the orbit is radially stable, we can approximate w S(i? Q ), discarding second-order and higher terms. Thus, 

$ eff (R, z) a $ e //Cft, 0) + 27tG£( j R )|*|, (28) 

and the corresponding approximate Hamiltonian, H = (Pj| + P^)/2 -I- $ e //, is also separable. Thus, since \z\ is 
continuous, we can solve the corresponding Hamilton- Jacobi equation assuming that the generating function has the 
form 

S(R, z, J) — Sr(R, J) + S z (z, J), (29) 

where J are the approximate action variables. In this way, we obtain by quadratures two independent integrals of 
motion near (i? o ,0): one for the i?-coordinate, Jr, and another for the z-coordinate, J z . In the next subsection we 
study with some detail this idea, focusing on the z-component, which is of special interest here. 



4.1. Shape of nearly equatorial orbits 
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Fig. 2. — Mass density profiles, S (10 10 kg m — 2 ) vs. R/a, for the nine spiral galaxies of figure[l]via relation 

Consider an orbit in the equatorial plane and the corresponding vertical perturbation. If the original orbit is circular, 
motion will occur on the torus given by the approximate action variables Jr, J z , which are obtained by solving the 
Hamilton- Jacobi equation with a separable generating function (|29[) : 



E = 



2 V dR 



2 V dz 



$ e// (i?,0) + 27rGS(i? o )|^| 



We have that 



is an integral of motion, where 
In consequence, 

with H z = E z . Thus S z is given by 



H s = ^ + u\z\ 



lu = 2ttGY 1 {R ). 



(30) 
(31) 



2 \ dz 



■ U)\Z\ = Ey 



S Z = V2~J y/E z -ui\z\dz, 



depending only on z. The corresponding action variable is (see iBinnev &: Tremain e (2008) for a discussion) J z = 
AS z /2ir, and in consequence, 

2\^2 f'^' 

J z = / y/E z - uj\z~\dz, 

n Jo 

where Z = z max is the maximum value of z along the orbit. From (1301) we have Z — E z /ui and, since z is positive 
along the range of integration, we can write J z — 4:y/2E z / Sttcj or, in terms of Z, 

V2 1/2 , 



Jz 



3tt 



. w l/2 Z 3/2^ 



(32) 
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Fig. 3. — Quadratic epicyclic frequency, k 2 (10 22 s 2 ) vs. R/a, for the sample of figures IT1 and l2l In all of these cases we obtain re 2 > 0, 
meaning that quasi-circular orbits, on the optical disk region, are stable under radial perturbations. 

The action-angle formalism described above for the z-part of the approximate Hamiltonian Q30p gives exactly the 
period (fTS]) and the amplitude (fl"9|) for the vertical oscillations, as expected. 

4.2. Adiabatic invariance of J z 

We now study the eff ects of adiabatic variations in the approximate potential of eq. ([30]) (see section 3.6 of 
iBinnev fc Tremaind ([2008)). If we consider a slow change in u> (see eqs. (1301) and (I3ip ). 



where s « 1 is a dimensionless quantity that varies slowly with time, adiabatic invariance of J z gives uj' x I 2 Z' 3 ! 2 — 
w 1 / 2 Z^l 2 , and we have Z' = Zs^ 1 / 3 , or equivalently, 

1/3 

(33) 



2T_ 
~Z 



Now consider R as a function of time, which induces an "effective" time variation in the parameter ui (defined by 
(j3"Tj) ) appearing in (|3"0")) : 

u){t)=2irGT.(R{t)). (34) 

Assuming that this effect is an adiabatic variation in the approximate potential of eq. (|30l) . corresponding to the 
perturbed equatorial orbit in the effective potential (l27l) . eq. (|33|) relates the maximum z-amplitudes of a given orbit 
at different values of R through the following expression: 



Z(R) 



Z(R) \Z(R) 



1/3 



(35) 



This relation determines the "envelope" of the orbit in the meridional plane. Therefore we expect that if the char- 
acteristic frequencies of the vertical perturbations are much higher than the frequencies of the i?-perturbations, then 
this approximation will be valid and the maximum z-amplitude of motion will be given by eq. (|35[) . This result can be 
compared with the "envelopes" of (numerically calculated) 3D orbits in a given potential $ e // near its minimum. Eq. 
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Fig. 4. — Kuzmin disk: Dimensionless effective potential in the equatorial plane for i/yaGm, = 0.2 (left side) and dimensionless surface 
density (right side). 



3) m ust also be compared with the corresponding relation for smooth potentials (eq. (3.279) of iBinnev fc Tremaind 
2008j)). 



4.3. Numerical Experiments 

In order to illustrate the applicability of the ideas sketched above, we perform simulations for the motion of test 
particles around mass distributions of the form ([8j) • We shall focus on the validity of relation (|35l) ; for this reason we 
present a number of situations in which the orbits integrated are far from being considered as "nearly equatorial" or 
"nearly circular" . 

At first, we present the results of numerically calculated orbits in two cases: (i) the Kuzmin disk (fig. [5J and (ii) 
the Kuzmin disk surrounded by a Plummer halo (fig. [7]). The APDP for the Kuzmin disk is given by 



Si. = 



Gm 



27TO 2 



1 R 2 



-3/2 



and the APDP for the Plummer halo is 



GM 



*p = — /rr 



3M 



Pp 



r 

4tt& 3 V 1 + 62 



b 2 

2\ -5/2 



(36) 
(37) 

(38) 
(39) 



where r 2 — R 2 + z 2 (|Binnev fc Tremafnel (|2008|) ). In the case (i) we choose l/\/aGm = 0.2 and aE/Gm = —0.45 for 
numerical integrations. For this value of angular momentum the effective potential has a minimum near R/a — 0.5 (see 
fig- El and for this value of energy we can obtain disk-crossing orbits in the region < R/a < 2. We find that there is 
a great number of orbits obeying the relation (|35|1 . even in regions far away from the critical point (i.e. corresponding 
to the radius of the circular orbit), where the approximation (|28j) presumably should not be valid. We remark that in 
all of these cases, the test particle passes through (or very near) the critical point, at least one time. Figure [5] shows 
two examples of this fact by plotting the motion in the meridional plane of two orbits along with the orbit's envelope 
computed from (1331) . Note that the prediction of ([33)) holds well beyond the vicinity of the thin disk, as it can be 
viewed in the right panel of fig. [5] 

In the case (ii), with the addition of the Plummer halo, we consider three situations: M/m = 0.01,0.5, 1, in order 
to account gradually the contribution of this component in the orbits' behavior. First we choose an extended halo 
(b/a = 2) and then a more concentrated one (b/a = 0.2), as it is shown in figures [6] and respectively. By maintaining 
the same values of energy and angular momentum as in the above case, we obtain the meridional-plane orbit for each 
ratio M/m. The results of the computation for b/a = 2 are shown in fig El from which we can see small deviations 
from the predictions of (|35p when M/m = 0.5 and M/m = 1 (for M/m = 0.01 the test particle describes a nearly 
equatorial orbit which is very well modeled by (|35p ). It is remarkable that such deviations are not very significant 
when the halo mass is of the order of the disk mass and the amplitude of vertical oscillations is comparable with the 
effective size of the disk. 

In contrast, the results of the computation for b/a — 0.2, which are shown in fig. [SI reveal significant deviations 
from the prediction of (1351) . This is due to the fact that, in all of the situations illustrated, the orbit does not pass 
through the critical point of the effective potential (see fig. [8]). However, note that for M/m = 0.01 the orbit passes 
(say) near the critical point and the deviation is significantly smaller than in the other two cases. 
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Fig. 5. — Motion in the meridional plane around a Kuzmin disk with aE/Gm = —0.45, t/\faGm = 0.2 and initial conditions z/a = 10 -15 , 
Pr = , R/a=0.2 (left), R/a=0.25 (right). The prediction of eq. H35I) is in red and the meridional-plane orbit is in black. The orbit in the 
left panel can be considered as a nearly equatorial orbit and its vertical amplitude is very well modeled by the red line. The same happens 
with the orbit in the right panel, although its vertical amplitude it is not near the equatorial plane. 
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Fig. 6. — Kuzmin disk + Plummer halo with b/a = 2: a) Effective potential in the equatorial plane, using e/VaGm = 0.2 and 
M/m = 0.01,0.5,1 (from top to bottom), b) Density of the Plummer halo for M/m = 0.01,0.5,1 (from bottom to top). The surface 
density of the Kuzmin disk is given in the right panel of fig. [4] 




Fig. 7. — Orbits in Kuzmin + Plummer potential, with aE/Gm = — 0.45 , £/V aGm = 0.2, initial conditions R/a = 0.2, z/a = 10 15 , 
Pr = and using the same parameters of Fig. [5] The prediction of eq. (1356 is in red for M/m = 0.01, 0.5, 1.0, from the left to the right. 
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Fig. 9. — Orbits in Kuzmin + Plummer potential, with aE/Gm = — 0,45 , tj\J aGra = 0.2, initial conditions R/a = 0.3, z/a = 10 15 , 
Pr = and using the same parameters of Fig. [8] The prediction of eq. (135 I t is in red for M/m = 0.01, 0.5, 1, from the left to the right. 

As a second example, we perform nu merical calculations with a mathematically more involved model: the second 
member of the generalized Kalnajs disks (Gonzalez & Rcina (2006)) immersed in a halo's spherical logarithmic poten- 
tial. The surface mass density of the thin disk is a monotonically decreasing function, £ = 5Af/(27ra 2 )(l — R 2 /a 2 ) 3 ^ 2 , 
where a is the radius of the disk (it is a disk with finite extension, as in subsection 13. 3|) and M its total mass. 
The corresponding gravitational potential can be cast in oblate spheroidal coordinates, £ = a~ 1 Re[y/ R 2 + (z — ia) 2 ], 
?7 = — a _1 Im[A/i? 2 + (z — ia) 2 ], through the relation 

+B(35t7 4 - 30?7 2 + 3)] , (40) 



<S>K2 = [cot- 1 £ + A{3 V 2 - 1) 



with 



^7j[( 3 e 2 + 1 ) cot ~ 1 £- 3 £]< ( 41a ) 



B = 



14 
3 
448 



(35f 4 + 30£ 2 + 3) cot" 1 £ - 35<f - — £ 



(41b) 



This potential leads to a Keplerian rotation curve, in contrast with the first member of the family (the well known 
Kalnajs disk), which describes a configuration rotating as a rigid body. On the other hand, the potential modeling the 
spherical halo is given by 

$ H = V 2 ln(i? 2 + z 2 + d 2 ) (42) 

where V and d are parameters to be determined according to the observation al data (for example, the galactic rotation 
curve). This component was used by Uohnston. Spiergel fc Hernquistl ([1995) to simulate the motion of dwarf galaxies 
around the Milky Way. Here we present the results of numerical experiments of the test-particle motion in the potential 
+ by considering two situations: (i) A case in which the disk component is dominant, where Vy/ a/GM = 0.1, 
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Fig. 10. — Model ® K2 + *H.Left panel: Dimensionless effective potential for ij\/aGM = 0.12, V y'a/GM = 0.1, d/a = 1.2 (lower 
curve) and lj\J aGM = 0.8, V \/a/GM = 1.2, d/a = 0.6 (upper curve). Right panel: Dimensionless Circular velocity for l/\J aGM = 0.12, 
V^/a/GM = 0.1, d/a = 1.2 (lower curve) and £/V 'aGM = 0.8, V y^a/GM = 1.2, d/a = 0.6 (upper curve). 



d/a = 1.2, for the halo, and £/VaGM = 0.12, Ea/GM = -1 and Ea/GM = -2, for the orbits; (ii) a case characterized 

by a dominant halo component, with V \J a/GM = 1.2, d/a = 0.6, and £/VaGM = 0.8, Ea/GM — —0.5 for orbits. 
The effective potential and rotation curve corresponding to both cases are shown in fig. 1101 

In the case (i) we obtain results similar to the reported ones for the Kuzmin-Plummer potential. Working with 
Ea/GM = —1, we find orbits which do not pass through the critical point of the effective potential and, in consequence, 
do not follow the prediction of eq. (|35|) (top pannels of fig. H"Tj) . In contrast, we find that for Ea/GM = —2, there is a 
great number of orbits whose vertical amplitude is described by (|35p with high precision. In the bottom panels of fig. 
fTTI we show two of them. 

In situation (ii) we note an interesting phenomenon. Here we find orbits which, in general, do not follow eq. (|35p 
but there is a number of them which pass near the red line where the particle reaches the maximum amplitude. This 
is the case illustrated in the top panels of fig. [12] where we present a chaotic orbit (left side) and a loop orbit (right 
side), which pass through the critical point of the effective potential (near R = 0.5) and also approach the red line. 
Note that the z-amplitudc of the oscillations is comparable with the i?-amplitude. The bottom panels of this figure 
show two examples of regular orbits which do not obey the prediction of |35|) . The chaoticity or regularity of the 
aforementioned orbits can be visualized in the surface of section of fig. 1131 

We point out that the accurateness of eq. (|35p for numerically integrated orbits in the above examples (specially 
when dealing with nearly equatorial orbits) can also be seen as an indirect verification of the stability criterion (|15p . 
since the derivation of eq. (|35p depends on the formalism introduced to obtain condition (|15p . 

We performed numerical integrations with the Runge-Kutta method of fourth order with variable time step. Con- 
servation of energy was checked with a precision characterized by a maximum relative error of about 10 -6 . All 
the computations were performed in the Laboratorio de Computagao Paralela Patricio Letelier at the Instituto de 
Matematica, Estati'stica e Computagao Cicntifica of UNICAMP. 

5. STABILITY CRITERIA IN MODIFIED THEORIES OF GRAVITY 

The results obtained for the vertical stability of equatorial circular orbits (as well as the predictions about the 
amplitude of nearly circular orbits) are valid for Newtonian gravity, but there is no guarantee that this also happens in 
other theories. However, most theories of modified gravity have the property that, in the Newtonian limit, test particles 
do follow a Hamiltonian flow determined by a modified potential. This makes it possible to formulate, in a similar 
fashion as in Newtonian gravity, a criterion for the Liapunov stability of circular orbits and an adequate description 
of the oscillations of nearly circular orbits (assuming adiabatic invarian ce of the approximate vertica l action). In 
this section we illustrate this f a ct by addressing two examples: MOND (Bekenstein & Milgrom (1984)) and RGGR 
(jRodrigues. Shapiro fc Letelier! (j2010f )). 

MO ND can be formulated as a potential theory, with a potential VE' satisfying the field equation (Bekenstein & Milgrom 

V ■ |)i(|V*|/o )V*] - 4ttGp, (43) 



where fj,(x) is the MOND interpolating function, 
density. For a thin disk of matter in the plane z 



a Q is a characteristic acceleration scale and p is the baryonic matter 
= we obtain from Gauss's theorem 



<9* 

W\ 



z=0 



2ttGY,{R) 

mOvvpik)' 



(44) 



which enables us to conclude that the condition for vertical stability of equatorial circular orbits is £ > 0, the same 
as in Newtonian gravity. Moreover, it can be verified that the condition for radial stability is k 2 > 0, where n 2 is 
calculated from the MOND potential 
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Fig. 11. — Orbits in the meridional plane for the model *K2 + $ff with 1/VaGM = 0.12, Vy/a/GM = 0.1, d/a = 1.2. Top panels: 
Ea/GM = -1, R/a = 0.2 (left) and R/a = 0.6 (right). Bottom panels Ea/GM = -2, R/a = 0.25 (left) and R/a = 0.35 (right). In all 
panels, z/a = 10 5 . 

On the other hand, we find that the adiabatic invariance of the approximate vertical action leads to the relation 

1/3 



Z{R) 



(45) 



which in terms of the baryonic surface mass density of the thin disk reads 

Z{R) = ( m Z(R) Y /3 

Z{R) \n k nR) J ' 

where ^ = /i(|W(JZ,0)|/a o ). 

For RGGR we have the modified potential (jRodrigues. Shapiro fc Letelierl ([2010)) 

c 2 

* = ®n + 2^G, 

where $jv is the Newtonian potential, Go is the "standard" gravitational constant and G = G{$>n) takes into account 
renormalization group corrections. The phenomenologically adopted form of G gives us (see Rodrigucs, Sha piro fc Letelierl 



(46) 
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Fig. 12.— Orbits in the meridional plane for the model $ X 2 + with t/VaGM = 0.8, V^/a/GM = 1.2, d/a = 0.6, Ea/GM = -0.5, 
and initial conditions R/a = 0.8, -R/a = 0.81, R/a = 0.78, -R/a = 0.76 (z/a = 10 -15 in all panels). Note the chaotic orbit which is bounded 
by the red line with a reasonably good precision. 



K3» 



v 2 

v OQ 

$7V 



JV, 



where is the asymptotic circular velocity and $jv -^0 as |x| — > oo. Therefore we have 



9$ 



2vrG 



2 = 



r 2 



AT 



(47) 



(48) 



and we see that the vertical stability condition for circular orbits is again £ > 0, and the radial stability condition is 
also n 2 > (computed with the potential $). Adiabatic invariance of J z for nearly circular orbits gives us (see eq. 

USD) 



z{R) 

Z(R') 



1 - 



S(i?') 



1 - 



E(iJ) 



1/3 



(49) 



*w(ii,0) / 

We see that modified theories of gravity predict deviations from the Newtonian behavior of nearly equatorial orbits. 
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Fig. 13. — Surface of section (z = 0, P z > 0) formed by the consequents of the four orbits in fig. 1121 

For orbits of astrophysical objects with a considerable radial amplitude these different predictions could, in principle, 
be compared with astronomical data. 

6. CONCLUSIONS 

We analyzed the stability of equatorial circular orbits in thin disks surrounded by smooth axisymmetric structures. 
The presence of the thin disk does not allow us to proceed in the same way as in the smooth case because of the delta- 
like singularity. In particular, the vertical stability criterion for smooth potentials, v 2 > 0, is not applicable anymore. 
We developed a consistent vertical stability criterion for circular orbits, which together with the (unchanged) radial 
stability criterion, ensures Liapunov stability of the corresponding circular orbit. Based on this new formalism, we find 
that nearly equatorial orbits have a third integral of motion, which is given by eq. (|35[) . This is supported by numerical 
simulations, which additionally reveal that orbits with great vertical amplitude can be described approximately by this 
integral (sec. 14. 3|) . It would be interesting to see if this dependence on the surface density is present in more realistic 
models, not described by a razor-thin disk, but incorporating a stellar distribution described by a thickened disk. 

The introduction o f the new stability criterion leads to the conc l usion t hat all of the thin disk models presented 
in iGonzalez fc Reinal (|2006D ; [Gonzalez. Plata-Plata fc Ramos-Carol ([2010ft ; IPedraza. Ramos-Caro fc Gonzalez! (|2008h 
are stable in a first approximation, contrary to the statements shown in such references. This fact urged us to obtain 
additional models in sec. 13.31 in order to show that Hunter's method, taking into account the stability criterion 
constructed here, is a powerful tool to model the maximum disk of a number of flat galaxies. Having tested the 
stability of orbits in this class of models, it would be interesting to carry out more conclusive stability analyses based 
on statistical mechanics considerations (i.e. perturbed solutions of Boltz mann equation, Toomre's criterion, etc) . We 
also have to point out that the stability analys es performed in references iRamos-Caro, Pedraza &: Letelieil ([20111 ) and 
iRamos-Caro. Lopes-Suspes fc G onzalez (|2008l) need to be corrected. 

We also briefly addressed the problem of the stability criterion in modified theories of gravity, such as MOND and 
RGGR, in sec. [5] (the same problem, in the realm of general relativity theory, is being studied and the results will be 
shown in a next paper) . We point out that whenever is possible establish a Hamiltonian formulation of the motion, it 
is also possible to perform an analysis similar to the presented here, in the Newtonian gravity realm. In particular, for 
the two examples studied here, we find that the stability criteria are also given by the relations £ > and k 2 > and 
that the assumption of adiabatic invariance leads to eqs. (|4"5j) and ([4"9")l , introducing deviations from the Newtonian 
relation (j35[) . This fact can be used as an additional test of these theories, once we have at disposal the required 
observational data. 
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APPENDIX 

PROOF OF THE VERTICAL STABILITY CONDITION 
The z- equation 

Consider a circular orbit of radius Rq in the plane z = 0, with z-component of angular momentum given by I. We 
define a vertical perturbation on this orbit as an instantaneous increase in the z-component of the velocity of the 
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particle (at time t = 0, say): 

v(0)^v{0) + v Qz z. (Al) 

This "shift" in velocity does not change the value of I, but increases the energy of the orbit by an amount ^(vq z ) 2 . 

We are interested in small vertical perturbations to the circular orbit. The term "small" will become clearer in the 
course of the proof. The questions we want to answer are the following: 

1. For sufficiently small vq z , what are the conditions that make the particle come back to the plane of the disk in a 
finite time? 

2. If the particle comes back, is the z-amplitude of motion compatible with the assumption of a small perturbation? 
In this case, is the R-variation negligible? 

3. After coming back and crossing the disk, what happens to the particle? Does it stay near the original circular 
orbit for t — > oo ? 

Does the particle come back to the disk? 
Without loss of generality, assume v$ z > 0. Then the particle goes to the region z > 0, and z = — <9$/<9z implies 

m = m- f ^(R(t'), z (t'))dt>, (A2) 

Jo oz 

while the equation for R gives 

R{t) - R(0) = - J* ^l(R(t'), z(t'))dt'. (A3) 

Since R(0) = and d^ e f f /8R(Rq, 0) = 0, the continuity of d$ e ff/dR implies that for sufficiently small z(t') we will 
have R(t) negligible, and then the assumption of constant R will be reasonable. This condition is satisfied for small 
enough vq z , at least for a short time. Let us assume for now R w R and study its consequences. We will come back 
to this issue later on. Equation (IA2|) reduces to 

m=vo,-J -^(Ro,z(t'))df. (A4) 

Since d$/d\z\ is continuous and d$>/d\z\(R , 0) > 0, take a > such that d$/d\z\(R , z) > for z £ [—a, a]. If we 
define 

<9$ 



a Ro = mm i —(R a ,z) \, (A5) 

zS[-a,a] lO\z\ J 

we have that ctR Q > by compacity. 

If z(t') < a for t' e [0,t], it follows that z(t') < (since z = -d$/dz = —d&/d\z\), and the vertical velocity of the 
particle tends to decrease. Eq. (|A4[) implies 

z(t) < v 0z - a Ro t. (A6) 

In particular, z(t) < vo z , which implies 

z(t) < v 0z t. (A7) 

Then, defining 

t c = —, (A8) 
we have that z(t) < a for t < t c . Integrating (IA6[) . we also have the inequality 

z(t)<v Qz t-^t 2 . (A9) 



The first instant of time in which the right-hand side of eq. (|A9p is equal to a is 

1 



t' c = 

aR 



''02 



V (vq z ) 2 - 2a Ro a 



(A10) 



if the discriminant is positive. We also have z(t) < a for t < t' c . 

Since a and a Ro depend only on the potential and not on the particular trajectory of the particle, we can fix these 
quantities by imposing the condition that a is sufficiently small to keep R small enough, in such a way that the 
i?-variation of the perturbed orbit can be neclected at least for a short time (we will see below that this is guaranteed 
for long times by the radial stability of the circular orbit). Doing this, we have from the definition of t c and from eq. 
[6j that there is a critical value wTjl such that if vq z < uqJ then there will be a time t £ [0, t c ] (depending on vq z ) 
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such that z(t) < 0. Since a n n > and z(t) < a by construction, the particle will have z(t) < for t > t, and since the 
upper estimate of z(t), eq. (|A6[) . decreases with t, i(t) will not approach zero for larger t, which implies the particle 
will hit the disk in a finite time. 

We can obtain from eq. (|A9j) an uppe r estimate for the total time interval the particle stays in the z > region. If 
we make the right hand side of eq. (|A9I) equal to zero (which implies z(t) < 0), we find 

At=^£. (All) 

We can also obtain an estimate for vq z : from eq. (| A9|) . a sufficient condition to have z(t) < a during the whole 
oscillation is 

v 0z t-^t 2 <a. (A12) 

This condition will be satisfied for all t if, and only if, vq z < ^/2aR a. Thus, we have that a lower limit for the critical 
value voz~ is 

WI > y/2 aRo a. (A13) 

That is, if vq z < ^J2aR a the particle will hit again the plane of the disk within the time interval At given by (|A11I) 
and the oscillation around the original circular orbit will have a vertical amplitude smaller than a. In this way, we 
have also answered the second question: the amplitude of the oscillation is less than a (given the above conditions), 
where a can be taken arbitrarily small. In fact, we can obtain an estimate for this amplitude: given vq z , it follows 
from e qs. (|A9[ ) and (IA11[) that the amplitude of the perturbation will be smaller than the value of the right-hand side 
of eq. (|A"9")) evaluated at At/ 2: 

W < (A14) 

The case vq z < is analogous, because of the ^-symmetry of the system. Thus we have answered the first question: 
For small enough |i>o z |, the condition for the particle come back to the plane of the disk in finite time is £(i?o) > 0. 

Disk-crossing and asymptotic behavior 

Assume that the i?-variation of the vertically perturbed orbit is negligible. This implies that the projection of the 
perturbed orbit on the z = plane is the original circular orbit, and then by conservation of the mechanical energy, 
the particle will hit the disk with a velocity with vertical component — vq z . 

The equation of motion for z implies a discontinuity in the particle's acceleration while crossing the disk, but its 
velocity is continuous. This implies that the particle will have a velocity — vq z just after crossing the disk, and since the 
system is symmetric with respect to the z — plane the particle will strike again the disk (now from the other side). 
Motion after crossing the disk will be analogous to the oscillation before crossing it, because of the z-symmetry of the 
potential and the nature of the new initial conditions. Thus, it follows that for sufficiently small \vq z \ the vertically 
perturbed orbit will r emain oscillating around the original cir cular orbit for t — > oo, with characteristic period given 
by 2 At (see eq. (lAlip ) and characteristic amplitude given by (|A14I) . 

In this sense, the vertical stability condition E(i?o) > derived in this section has, for thin disks represented by 
density distributions of the form (151), the same status as the cond i tion v 2 > for smooth axisymmetric potentials 
(jBinney fc Tremaind (|2008f ) ; IGonzalez. Plata-Plata fc Ramos-Carol (|2010l ) ) . Both conditions assume small perturba- 
tions and neglect the i?-variation. We now analyze the effect of this variation and obtain a stability condition under 
general small perturbations. 

Effects on the R- coordinate of the trajectory 

In order to be able to neglect the R- variation of the perturbed orbit, a sufficient condition is that the original circular 
orbit is stable under small radial perturbations, as we shall see in the following. This translates into the inequality 

k 2 {R ) > 0, (A15) 

where k, 2 (Rq) = ^-^t^(Rq,0) is the quadratic epicyclic frequency of the perturbation (jBinnev fc Tremainel (|2008[ k 
IGonzalez. Plata-Plata fc Ramos-Carol (120101 )). Indeed, assuming S(i?o) > and k 2 (Rq) > (with Rq not on the 
border of the thin surface distribution if this distribution is finite in size) , we have 

e// -(i? ,0)>0, (A16) 



OR 2 
d\z\ 



CRo,0)>0, (A17) 



which imply (Ro,0) is a (strict) local minimum of the potential & e ff(R,z): there is a neighborhood of (i?o,0) in 
which $ e ff(R,z) > <& e //(i?o, 0) if (R, z) ^ (Rq,0). Thus, for small enough vq z (corresponding to small enough 
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(E — QefflRo, 0))), the orbit will oscillate around (Rq, 0) in the meridional plane with an amplitude that can be made 
arbitrarily small. 

It also follows from conditions (|A16I) and (|A17[) that, if we neglect changes in I due to small radial perturbations, 
the circular orbit will be Liapunov stable under small perturbations in any direction of the meri dional p l ane. In fact, 
Liapunov stability of (Rq, 0) depends only on the continuity of & e ff and not on its smoothness (see lArnoldl ()1978[ ). chap. 
5, p. 99). This result is the generalization to thin dis ks represented by density di s tributions of the form Q of the general 
stabilit y condition v 2 > 0, k 2 > for circu lar orbits (jBinnev fc Tremaind (|2008f k [Gonzalez. Plata-Plata fc Ramos-Carol 
(p?JTfl : lRamos-Caro. Pedraza fc Leteherl (poll ). 

Finally we note that, while the condition S(i?o) > depends only on the thin disk, condition k 2 (Rq) > depends 
also on the smooth 3D distribution, in such a way that this spatial distribution can affect the radial stability of the 
circular orbit and, as a consequence, the behavior of vertically perturbed orbits for large enough time. 
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